
function G = grad_cs_moments(gamma,B,X0,C,xi0,Z,KD,n,optInner,cmax,rfd)

N0 = n(1,1); N1 = n(1,2); N2 = n(1,3); N = N0 + N1 + N2;

% gradient of campaign_stage moment conditions (by central finite differences)
G = zeros(size(Z,2),length(gamma)+length(B));
% w.r.t. party payoff parameters
for k = 1:length(gamma)
    tic
    df = 10^round(log10(rfd * mean(abs(gamma)))); % step size
    gamma_df = gamma;
    gamma_df(k) = gamma_df(k) - df;
    parfor t = 1:N0
        tt = [t,N+t,2*N+t,2*N+N0+N2+t,3*N+N0+t];
        br0(t,:) = BRall(gamma_df,B,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,0); %#ok<*PFBNS>
    end
    parfor t = 1:N1
        tt = [N0+t,N+N0+t,1,2*N+N0+N2+N0+t,3*N+N0+N0+t];
        br1(t,:) = BRall(gamma_df,B,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,1);
    end
    parfor t = 1:N2
        tt = [N0+N1+t,N+N0+N1+t,2*N+N0+t,1,3*N+N0+N0+N1+t];
        br2(t,:) = BRall(gamma_df,B,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,2);
    end
    br = [br0;br1;br2];
    br = [br(1:N,1);br(1:N,2);br([1:N0,N0+N1+(1:N2)],3);br(1:(N0+N1),4);br(1:N,5)];
    g1 = Z' * (C - br);
    gamma_df = gamma;
    gamma_df(k) = gamma_df(k) + df;
    parfor t = 1:N0
        tt = [t,N+t,2*N+t,2*N+N0+N2+t,3*N+N0+t];
        br0(t,:) = BRall(gamma_df,B,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,0);
    end
    parfor t = 1:N1
        tt = [N0+t,N+N0+t,1,2*N+N0+N2+N0+t,3*N+N0+N0+t];
        br1(t,:) = BRall(gamma_df,B,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,1);
    end
    parfor t = 1:N2
        tt = [N0+N1+t,N+N0+N1+t,2*N+N0+t,1,3*N+N0+N0+N1+t];
        br2(t,:) = BRall(gamma_df,B,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,2);
    end
    br = [br0;br1;br2];
    br = [br(1:N,1);br(1:N,2);br([1:N0,N0+N1+(1:N2)],3);br(1:(N0+N1),4);br(1:N,5)];
    g2 = Z' * (C - br);    
    G(:,k) = (g2 - g1) / (2 * df);    
    toc
    disp('G gamma')
    disp(k)
end
% w.r.t. demand menu-party fixed effects
% independent candidates
tic
kk = [1,3+KD+1,2*(3+KD)+1,2*(3+KD)+2+KD+1,2*(3+KD)+2*(2+KD)+1];
df = 10^round(log10(rfd * mean(abs(B(kk)))));
parfor t = 1:N1
    tt = [N0+t,N+N0+t,1,2*N+N0+N2+N0+t,3*N+N0+N0+t];
    br1(t,:) = BRall(gamma,B,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,1);
end
parfor t = 1:N2
    tt = [N0+N1+t,N+N0+N1+t,2*N+N0+t,1,3*N+N0+N0+N1+t];
    br2(t,:) = BRall(gamma,B,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,2);
end
for k = kk
    B_df = B;
    B_df(k) = B_df(k) - df;
    parfor t = 1:N0
        tt = [t,N+t,2*N+t,2*N+N0+N2+t,3*N+N0+t];
        br0(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,0);
    end
    br = [br0;br1;br2];
    br = [br(1:N,1);br(1:N,2);br([1:N0,N0+N1+(1:N2)],3);br(1:(N0+N1),4);br(1:N,5)];
    g1 = Z' * (C - br);
    B_df = B;
    B_df(k) = B_df(k) + df;
    parfor t = 1:N0
        tt = [t,N+t,2*N+t,2*N+N0+N2+t,3*N+N0+t];
        br0(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,0);
    end
    br = [br0;br1;br2];
    br = [br(1:N,1);br(1:N,2);br([1:N0,N0+N1+(1:N2)],3);br(1:(N0+N1),4);br(1:N,5)];
    g2 = Z' * (C - br);    
    G(:,length(gamma)+k) = (g2 - g1) / (2 * df);
    toc
    disp('G beta')
    disp(length(gamma)+k)
end
% PRI coalition candidate
tic
kk = [2,3+KD+2,2*(3+KD)+2+KD+2,2*(3+KD)+2*(2+KD)+2];
df = 10^round(log10(rfd * mean(abs(B(kk)))));
parfor t = 1:N0
    tt = [t,N+t,2*N+t,2*N+N0+N2+t,3*N+N0+t];
    br0(t,:) = BRall(gamma,B,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,0); %#ok<*PFBNS>
end
for k = kk
    B_df = B;
    B_df(k) = B_df(k) - df;
    parfor t = 1:N1
        tt = [N0+t,N+N0+t,1,2*N+N0+N2+N0+t,3*N+N0+N0+t];
        br1(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,1);
    end
    br = [br0;br1;br2];
    br = [br(1:N,1);br(1:N,2);br([1:N0,N0+N1+(1:N2)],3);br(1:(N0+N1),4);br(1:N,5)];
    g1 = Z' * (C - br);
    B_df = B;
    B_df(k) = B_df(k) + df;
    parfor t = 1:N1
        tt = [N0+t,N+N0+t,1,2*N+N0+N2+N0+t,3*N+N0+N0+t];
        br1(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,1);
    end
    br = [br0;br1;br2];
    br = [br(1:N,1);br(1:N,2);br([1:N0,N0+N1+(1:N2)],3);br(1:(N0+N1),4);br(1:N,5)];
    g2 = Z'*(C-br);    
    G(:,length(gamma)+k) = (g2 - g1) / (2 * df);
    toc
    disp('G beta')
    disp(length(gamma)+k)
end
% PVEM coalition candidate
tic
kk = [3,3+KD+3,2*(3+KD)+2,2*(3+KD)+2*(2+KD)+3];
df = 10^round(log10(rfd * mean(abs(B(kk)))));
parfor t = 1:N1
    tt = [N0+t,N+N0+t,1,2*N+N0+N2+N0+t,3*N+N0+N0+t];
    br1(t,:) = BRall(gamma,B,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,1);
end
for k = kk
    B_df = B;
    B_df(k) = B_df(k) - df;
    parfor t = 1:N2
        tt = [N0+N1+t,N+N0+N1+t,2*N+N0+t,1,3*N+N0+N0+N1+t];
        br2(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,2);
    end
    br = [br0;br1;br2];
    br = [br(1:N,1);br(1:N,2);br([1:N0,N0+N1+(1:N2)],3);br(1:(N0+N1),4);br(1:N,5)];
    g1 = Z' * (C - br);
    B_df = B;
    B_df(k) = B_df(k) + df;
    parfor t = 1:N2
        tt = [N0+N1+t,N+N0+N1+t,2*N+N0+t,1,3*N+N0+N0+N1+t];
        br2(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,2);
    end
    br = [br0;br1;br2];
    br = [br(1:N,1);br(1:N,2);br([1:N0,N0+N1+(1:N2)],3);br(1:(N0+N1),4);br(1:N,5)];
    g2 = Z' * (C - br);    
    G(:,length(gamma)+k) = (g2 - g1) / (2 * df);
    toc
    disp('G beta')
    disp(length(gamma)+k)
end
% w.r.t. demand demographics
for k = 1:KD
    % MP
    tic
    df = 10^round(log10(rfd * mean(abs(B([3,3+KD+3,2*(3+KD)+2,2*(3+KD)+2+KD+2,2*(3+KD)+2*(2+KD)+3] + k)))));
    B_df = B;
    B_df(3+k) = B_df(3+k) - df;
    parfor t = 1:N0
        tt = [t,N+t,2*N+t,2*N+N0+N2+t,3*N+N0+t];
        br0(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,0);
    end
    parfor t = 1:N1
        tt = [N0+t,N+N0+t,1,2*N+N0+N2+N0+t,3*N+N0+N0+t];
        br1(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,1);
    end
    parfor t = 1:N2
        tt = [N0+N1+t,N+N0+N1+t,2*N+N0+t,1,3*N+N0+N0+N1+t];
        br2(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,2);
    end
    br = [br0;br1;br2];
    br = [br(1:N,1);br(1:N,2);br([1:N0,N0+N1+(1:N2)],3);br(1:(N0+N1),4);br(1:N,5)];
    g1 = Z' * (C - br);
    B_df = B;
    B_df(3+k) = B_df(3+k) + df;
    parfor t = 1:N0
        tt = [t,N+t,2*N+t,2*N+N0+N2+t,3*N+N0+t];
        br0(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,0);
    end
    parfor t = 1:N1
        tt = [N0+t,N+N0+t,1,2*N+N0+N2+N0+t,3*N+N0+N0+t];
        br1(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,1);
    end
    parfor t = 1:N2
        tt = [N0+N1+t,N+N0+N1+t,2*N+N0+t,1,3*N+N0+N0+N1+t];
        br2(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,2);
    end
    br = [br0;br1;br2];
    br = [br(1:N,1);br(1:N,2);br([1:N0,N0+N1+(1:N2)],3);br(1:(N0+N1),4);br(1:N,5)];
    g2 = Z' * (C - br);    
    G(:,length(gamma)+3+k) = (g2 - g1) / (2 * df);    
    toc
    disp('G beta')
    disp(length(gamma)+3+k)
    % NA
    tic
    B_df = B;
    B_df(3+KD+3+k) = B_df(3+KD+3+k) - df;
    parfor t = 1:N0
        tt = [t,N+t,2*N+t,2*N+N0+N2+t,3*N+N0+t];
        br0(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,0);
    end
    parfor t = 1:N1
        tt = [N0+t,N+N0+t,1,2*N+N0+N2+N0+t,3*N+N0+N0+t];
        br1(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,1);
    end
    parfor t = 1:N2
        tt = [N0+N1+t,N+N0+N1+t,2*N+N0+t,1,3*N+N0+N0+N1+t];
        br2(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,2);
    end
    br = [br0;br1;br2];
    br = [br(1:N,1);br(1:N,2);br([1:N0,N0+N1+(1:N2)],3);br(1:(N0+N1),4);br(1:N,5)];
    g1 = Z' * (C - br);
    B_df = B;
    B_df(3+KD+3+k) = B_df(3+KD+3+k) + df;
    parfor t = 1:N0
        tt = [t,N+t,2*N+t,2*N+N0+N2+t,3*N+N0+t];
        br0(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,0);
    end
    parfor t = 1:N1
        tt = [N0+t,N+N0+t,1,2*N+N0+N2+N0+t,3*N+N0+N0+t];
        br1(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,1);
    end
    parfor t = 1:N2
        tt = [N0+N1+t,N+N0+N1+t,2*N+N0+t,1,3*N+N0+N0+N1+t];
        br2(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,2);
    end
    br = [br0;br1;br2];
    br = [br(1:N,1);br(1:N,2);br([1:N0,N0+N1+(1:N2)],3);br(1:(N0+N1),4);br(1:N,5)];
    g2 = Z' * (C - br);    
    G(:,length(gamma)+3+KD+3+k) = (g2 - g1) / (2 * df);    
    toc
    disp('G beta')
    disp(length(gamma)+3+KD+3+k)
    % PVEM
    tic
    B_df = B;
    B_df(2*(3+KD)+2+k) = B_df(2*(3+KD)+2+k) - df;
    parfor t = 1:N0
        tt = [t,N+t,2*N+t,2*N+N0+N2+t,3*N+N0+t];
        br0(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,0);
    end
    parfor t = 1:N1
        tt = [N0+t,N+N0+t,1,2*N+N0+N2+N0+t,3*N+N0+N0+t];
        br1(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,1);
    end
    parfor t = 1:N2
        tt = [N0+N1+t,N+N0+N1+t,2*N+N0+t,1,3*N+N0+N0+N1+t];
        br2(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,2);
    end
    br = [br0;br1;br2];
    br = [br(1:N,1);br(1:N,2);br([1:N0,N0+N1+(1:N2)],3);br(1:(N0+N1),4);br(1:N,5)];
    g1 = Z' * (C - br);
    B_df = B;
    B_df(2*(3+KD)+2+k) = B_df(2*(3+KD)+2+k) + df;
    parfor t = 1:N0
        tt = [t,N+t,2*N+t,2*N+N0+N2+t,3*N+N0+t];
        br0(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,0);
    end
    parfor t = 1:N1
        tt = [N0+t,N+N0+t,1,2*N+N0+N2+N0+t,3*N+N0+N0+t];
        br1(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,1);
    end
    parfor t = 1:N2
        tt = [N0+N1+t,N+N0+N1+t,2*N+N0+t,1,3*N+N0+N0+N1+t];
        br2(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,2);
    end
    br = [br0;br1;br2];
    br = [br(1:N,1);br(1:N,2);br([1:N0,N0+N1+(1:N2)],3);br(1:(N0+N1),4);br(1:N,5)];
    g2 = Z' * (C - br);    
    G(:,length(gamma)+2*(3+KD)+2+k) = (g2 - g1) / (2 * df);    
    toc
    disp('G beta')
    disp(length(gamma)+2*(3+KD)+2+k)
    % PRI
    tic
    B_df = B;
    B_df(2*(3+KD)+2+KD+2+k) = B_df(2*(3+KD)+2+KD+2+k) - df;
    parfor t = 1:N0
        tt = [t,N+t,2*N+t,2*N+N0+N2+t,3*N+N0+t];
        br0(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,0);
    end
    parfor t = 1:N1
        tt = [N0+t,N+N0+t,1,2*N+N0+N2+N0+t,3*N+N0+N0+t];
        br1(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,1);
    end
    parfor t = 1:N2
        tt = [N0+N1+t,N+N0+N1+t,2*N+N0+t,1,3*N+N0+N0+N1+t];
        br2(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,2);
    end
    br = [br0;br1;br2];
    br = [br(1:N,1);br(1:N,2);br([1:N0,N0+N1+(1:N2)],3);br(1:(N0+N1),4);br(1:N,5)];
    g1 = Z' * (C - br);
    B_df = B;
    B_df(2*(3+KD)+2+KD+2+k) = B_df(2*(3+KD)+2+KD+2+k) + df;
    parfor t = 1:N0
        tt = [t,N+t,2*N+t,2*N+N0+N2+t,3*N+N0+t];
        br0(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,0);
    end
    parfor t = 1:N1
        tt = [N0+t,N+N0+t,1,2*N+N0+N2+N0+t,3*N+N0+N0+t];
        br1(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,1);
    end
    parfor t = 1:N2
        tt = [N0+N1+t,N+N0+N1+t,2*N+N0+t,1,3*N+N0+N0+N1+t];
        br2(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,2);
    end
    br = [br0;br1;br2];
    br = [br(1:N,1);br(1:N,2);br([1:N0,N0+N1+(1:N2)],3);br(1:(N0+N1),4);br(1:N,5)];
    g2 = Z' * (C - br);    
    G(:,length(gamma)+2*(3+KD)+2+KD+2+k) = (g2 - g1) / (2 * df);    
    toc
    disp('G beta')
    disp(length(gamma)+2*(3+KD)+2+KD+2+k)
    % PAN
    tic
    B_df = B;
    B_df(2*(3+KD)+2*(2+KD)+3+k) = B_df(2*(3+KD)+2*(2+KD)+3+k) - df;
    parfor t = 1:N0
        tt = [t,N+t,2*N+t,2*N+N0+N2+t,3*N+N0+t];
        br0(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,0);
    end
    parfor t = 1:N1
        tt = [N0+t,N+N0+t,1,2*N+N0+N2+N0+t,3*N+N0+N0+t];
        br1(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,1);
    end
    parfor t = 1:N2
        tt = [N0+N1+t,N+N0+N1+t,2*N+N0+t,1,3*N+N0+N0+N1+t];
        br2(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,2);
    end
    br = [br0;br1;br2];
    br = [br(1:N,1);br(1:N,2);br([1:N0,N0+N1+(1:N2)],3);br(1:(N0+N1),4);br(1:N,5)];
    g1 = Z' * (C - br);
    B_df = B;
    B_df(2*(3+KD)+2*(2+KD)+3+k) = B_df(2*(3+KD)+2*(2+KD)+3+k) + df;
    parfor t = 1:N0
        tt = [t,N+t,2*N+t,2*N+N0+N2+t,3*N+N0+t];
        br0(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,0);
    end
    parfor t = 1:N1
        tt = [N0+t,N+N0+t,1,2*N+N0+N2+N0+t,3*N+N0+N0+t];
        br1(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,1);
    end
    parfor t = 1:N2
        tt = [N0+N1+t,N+N0+N1+t,2*N+N0+t,1,3*N+N0+N0+N1+t];
        br2(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,2);
    end
    br = [br0;br1;br2];
    br = [br(1:N,1);br(1:N,2);br([1:N0,N0+N1+(1:N2)],3);br(1:(N0+N1),4);br(1:N,5)];
    g2 = Z' * (C - br);    
    G(:,length(gamma)+2*(3+KD)+2*(2+KD)+3+k) = (g2 - g1) / (2 * df);    
    toc
    disp('G beta')
    disp(length(gamma)+2*(3+KD)+2*(2+KD)+3+k)
end
% w.r.t. demand lagged vote share and campaign spending
for k = 1:3
    tic
    df = 10^round(log10(rfd * abs(B(3*(3+KD)+2*(2+KD) + k))));
    B_df = B;
    B_df(3*(3+KD)+2*(2+KD)+k) = B_df(3*(3+KD)+2*(2+KD)+k) - df;
    parfor t = 1:N0
        tt = [t,N+t,2*N+t,2*N+N0+N2+t,3*N+N0+t];
        br0(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,0);
    end
    parfor t = 1:N1
        tt = [N0+t,N+N0+t,1,2*N+N0+N2+N0+t,3*N+N0+N0+t];
        br1(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,1);
    end
    parfor t = 1:N2
        tt = [N0+N1+t,N+N0+N1+t,2*N+N0+t,1,3*N+N0+N0+N1+t];
        br2(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,2);
    end
    br = [br0;br1;br2];
    br = [br(1:N,1);br(1:N,2);br([1:N0,N0+N1+(1:N2)],3);br(1:(N0+N1),4);br(1:N,5)];
    g1 = Z' * (C - br);
    B_df = B;
    B_df(3*(3+KD)+2*(2+KD)+k) = B_df(3*(3+KD)+2*(2+KD)+k) + df;
    parfor t = 1:N0
        tt = [t,N+t,2*N+t,2*N+N0+N2+t,3*N+N0+t];
        br0(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,0);
    end
    parfor t = 1:N1
        tt = [N0+t,N+N0+t,1,2*N+N0+N2+N0+t,3*N+N0+N0+t];
        br1(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,1);
    end
    parfor t = 1:N2
        tt = [N0+N1+t,N+N0+N1+t,2*N+N0+t,1,3*N+N0+N0+N1+t];
        br2(t,:) = BRall(gamma,B_df,X0(tt,:),C(tt,:),xi0(tt,:),optInner,cmax,2);
    end
    br = [br0;br1;br2];
    br = [br(1:N,1);br(1:N,2);br([1:N0,N0+N1+(1:N2)],3);br(1:(N0+N1),4);br(1:N,5)];
    g2 = Z' * (C - br);    
    G(:,length(gamma)+3*(3+KD)+2*(2+KD)+k) = (g2 - g1) / (2 * df);
    toc
    disp('G beta')
    disp(length(gamma)+3*(3+KD)+2*(2+KD)+k)
end



